2.1.1 RasterLayer

# Chunk 1
# #1
e <- ext(c("xmin" = 34, "xmax" = 36, "ymin" = -6, "ymax" = -4))
#
# #2
r <- rast(x = e, res = 0.25, crs = crs(districts))
#
# #3            
set.seed(1)  
values(r) <- sample(1:100, size = ncell(r), replace = TRUE)  # 3
# r[] <- sample(1:100, size = ncell(r), replace = TRUE) 
# r <- setValues(x = r, values = sample(1:100, size = ncell(r), replace = TRUE))

par(mar = c(0, 0, 0, 4))
plot(st_geometry(districts), col = "grey", reset = FALSE)
plot(r, add = TRUE, ext = ext(districts))

We just used several functions from the terra package to create a random SpatRaster named r that has a 1/4 degree resolution and covers an area of 2 X 2 degrees in centural Tanzania. This particular raster is a temporary one that lives in memory.

Let’s walk through the labelled code. In # 1, we use terra’s ext (extent) function to define the boundaries of the raster, and then in # 2 use the rast function to create a raster from the resulting SpatExtent object e, assigning a CRS using the “crs” argument, which in turn uses terra’s crs to extract the crs from districts (#3). crs is similar to sf::st_crs, but outputs a different class of object that can’t be used by terra. The terra function can create a raster from many different types of input objects (passed to argument “x”), per ?rast:

filename (character), missing, SpatRaster, SpatRasterDataset, SpatExtent, SpatVector, matrix, array, list of SpatRaster objects. For other types it will be attempted to create a SpatRaster via (‘as(x, “SpatRaster”)’

Line # 2 creates an empty raster r with no cell values, so in # 3 we assign some randomly selected values into r. Note the method of assignment, using the values function; there are two other lines commented out below # 3 that show different ways of doing the same job.

The plot of r over districts uses the plot method defined for raster* objects. Note that it automatically creates a continuous legend, and also note that terra::plot can work with sf::plot. The plot is not very nice.

Let’s look now at the structure of the object r.

# Chunk 2
r
#> class       : SpatRaster 
#> size        : 8, 8, 1  (nrow, ncol, nlyr)
#> resolution  : 0.25, 0.25  (x, y)
#> extent      : 34, 36, -6, -4  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> name        : lyr.1 
#> min value   :     1 
#> max value   :   100
class(r)
#> [1] "SpatRaster"
#> attr(,"package")
#> [1] "terra"
typeof(r)
#> [1] "S4"
slotNames(r)
#> [1] "pntr"
r@pntr
#> C++ object <0x7664856c0> of class 'SpatRaster' <0x75a274200>
names(r@pntr)
#>  [1] "hasWindow"       "is_multidim"     "subset"          "range_max"       "hasCategories"  
#>  [6] "finalize"        "extent"          "timestep"        "initialize"      "hasUnit"        
#> [11] "names"           ".module"         "filenames"       "rgbtype"         "rgb"            
#> [16] "dataType"        "setValueType"    ".cppclass"       "hasRange"        "ncol"           
#> [21] "messages"        "hasColors"       "timezone"        "depth"           "origin"         
#> [26] "nsrc"            "hasDepth"        ".self"           "inMemory"        "get_sourcenames"
#> [31] "has_error"       "time"            "has_warning"     "hasValues"       "isMD"           
#> [36] "nlyr"            "range_min"       "valueType"       "get_crs"         "hasTime"        
#> [41] ".pointer"        "nrow"            "setValues"       ".refClassDef"    "res"            
#> [46] "units"
res(r)
#> [1] 0.25 0.25

r is an S4 object that has only slot pntr, which is actually from C++ (terra is written in C++ and creates a C++ object). There are a number of slots in there, which we can see by running names(r@pntr), but the terra packages provides methods for accessing those (e.g. ext() to get extent; sources() to get the raster’s filename, if written to disk; values() gets access to the data in the raster).

2.1.2 From 2- to 3-D

We have just seen how to create a SpatLayer and learned a bit about the structure of this kind of object, which is two-dimensional. We are now going to learn about three-dimensional rasters. Let’s create some new data.

# Chunk 3
r2 <- r > 50
r3 <- r
set.seed(1)
values(r3) <- rnorm(n = ncell(r3), mean = 10, sd = 2)
l <- list(r, r2, r3)

We used r to create two new rasters, r2 and r3. r2 was made by using a logical operator (>) to find the locations where r’s values exceed 50, creating a binary SpatRaster where 1 indicates the matching pixels, and 0 those that don’t. r3 was made by using r as a template, then overwriting the values with numbers generated randomly from a normal distribution (rnorm). These were then combined into list l.

# Chunk 4
s <- rast(l)
# s <- c(r, r2, r3)  # also works
names(s) <- c("r", "r2", "r3")
s
#> class       : SpatRaster 
#> size        : 8, 8, 3  (nrow, ncol, nlyr)
#> resolution  : 0.25, 0.25  (x, y)
#> extent      : 34, 36, -6, -4  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> names       :   r,    r2,       r3 
#> min values  :   1, FALSE,  5.57060 
#> max values  : 100,  TRUE, 14.80324
plot(s, nr = 1)

In the preceding code blocks, we created a multi-layer SpatRast from l, which is a series of rasters that have the same extent and resolution, which are “stacked” on top of one another in the order that they are given in the input list. The stacked layers can come from rasters held in memory, or from any number of files stored in different areas on disk, or from a multi-layer raster held in a single file on disk.

Applying plot to a stack or brick results in the automatic plotting of each layer into its own sub-window, with coordinates along the plot axes. You can specify numbers of rows and columns in your plotting device using the “nc” and “nr” arguments to terra::plot.

Here’s a more informative (by putting it over a map of Tanzania) way of plotting the layers in s. Note the use of ext in the plot, which we use to make the raster plot extent by the same as districts, to force the plot outside of district extent (this doesn’t work for the logical raster, however):

# Chunk 5
par(mfrow = c(1, 3), mar = c(0, 0, 0, 4))
for(i in 1:nlyr(s)) {
  districts %>% st_geometry %>% plot(col = "grey")
  plot(s[[i]], add = TRUE, ext = ext(districts))
}

2.3.1 Vector to raster

We have several vector datasets that come with sdadata which we can rasterize, starting with the district boundaries.

# Chunk 7
# #1
tzar <- rast(x = ext(districts), crs = crs(districts), res = 0.1)
values(tzar) <- 1:ncell(tzar)
#
# #2
districts <- districts %>% mutate(ID = 1:nrow(.))
distsr <- districts %>% 
  rasterize(x = ., y = tzar, field = "ID") %>% 
  print()
#> class       : SpatRaster 
#> size        : 108, 109, 1  (nrow, ncol, nlyr)
#> resolution  : 0.09958904, 0.09980746  (x, y)
#> extent      : 29.58953, 40.44473, -11.76235, -0.983143  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> name        :  ID 
#> min value   :   1 
#> max value   : 170

par(mar = c(0, 0, 0, 4))
plot(distsr, axes = FALSE)

In #1, we took an initial step to define a raster (tzar) that has the properties of resolution (0.1 decimal degrees), CRS, and extent that we want our rasterized vector to have. We set the extent of this raster to that of districts, using extent to get the bounding box coordinates (extent retrieves the same parameters as sf::st_bbox, but the output is in a different format).

In #2, we use rasterize to (as the name says) rasterize districts. The “y” argument is where we feed in tzar, the raster object that is the “target” for rasterizing districts. The “field” argument supplies the column names of the values that we want rasterized. In this case, we have to first create an “ID” variable for districts, in order to give an integer for each district name, as the district names (a character variable) cannot be written to the raster. Notice that we have constructed this as a pipeline (with print as the last step, to show the raster metadata). The commented out code below it shows how it can be done with more conventional syntax.

Our plot removes the coordinate-labelled axes and box that otherwise drawn around raster plots by default.

Next we rasterize the species dataset, which requires a little more prep to be meaningful:

# Chunk 8
# #1
tzar2 <- rast(x = ext(districts), crs = crs(districts), res = 0.25)
values(tzar2) <- 1:ncell(tzar2)

# #2
speciesr <- species %>% 
  dplyr::select(x, y) %>% 
  mutate(count = 1) %>% 
  st_as_sf(coords = c("x", "y"), crs = 4326) %>% 
  rasterize(x = ., y = tzar2, field = "count", fun = sum) %>% 
  print()
#> class       : SpatRaster 
#> size        : 43, 43, 1  (nrow, ncol, nlyr)
#> resolution  : 0.2524466, 0.2506792  (x, y)
#> extent      : 29.58953, 40.44473, -11.76235, -0.983143  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> name        :  sum 
#> min value   :    1 
#> max value   : 1216

# #3
par(mar = c(0, 0, 1, 4))
districts %>% 
  st_union %>% 
  plot(col = "grey", border = "grey", 
       main = expression(paste("N species per 0.25", degree, " cell")))
#> although coordinates are longitude/latitude, st_union assumes that they are planar
plot(speciesr, add = TRUE, axes = FALSE)

#> Error in match.names(clabs, names(xi)) : 
#>   names do not match previous names

In #1, we create a coarser scale (0.25 degree) version of tzar (tzar2), because we want our final raster to count how many farmers responding to our survey are found in each grid cell. The resolution of 0.1 degrees used for tzar2 is a bit too fine to convey the information nicely in a plot.

In #2, we do the rasterization as part of a pipeline. The first two lines prepare the species dataset. We add a new count variable (which assigns a value of 1 to each species record) before converting species to sf in the third line. We then rasterize the count variable, using the sum function to aggregate the number of species per 0.25\(\circ\) grid cell.

The resulting plot (#3) shows that most cells have less than 20 species. We add an extra plot decoration step, using expression with paste to add a superscript degree symbol to our plot title.

You can also rasterize line features, much as we did for points and polygons, which is shown below, but not run. The logic for not running it is that, with the previous raster packages, rasterizing lines was exceedingly slow, but it is actually quite fast with terra, so now we don’t run it because the rasterized lines aren’t that interested. However, there is code below that shows how one could so that with the roads dataset (subset to roads greater than 100 km long):

# Chunk 9
roadsr <- roads %>% 
  filter(as.numeric(st_length(.)) > 50000) %>% 
  mutate(ID = 1:nrow(.)) %>% 
  st_transform(crs = 4326) %>% 
  rasterize(., tzar, field = "ID")

2.3.2 Raster to vector

terra gives us functions that allow us to transform rasters to vectors, in this case terra’s native vector class, SpatVector, which we can further coerce to sf

# Chunk 10
# #1
dists_fromr <- as.polygons(x = distsr, dissolve = TRUE) %>% 
  print()
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 169, 1  (geometries, attributes)
#>  extent      : 29.58953, 40.44473, -11.76235, -0.983143  (xmin, xmax, ymin, ymax)
#>  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#>  names       :    ID
#>  type        : <int>
#>  values      :     1
#>                    2
#>                    3

dists_fromr <- st_as_sf(dists_fromr) %>% print()
#> Simple feature collection with 169 features and 1 field
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: 29.58953 ymin: -11.76235 xmax: 40.44473 ymax: -0.983143
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>    ID                       geometry
#> 1   1 MULTIPOLYGON (((36.66035 -3...
#> 2   2 POLYGON ((36.56076 -3.37852...
#> 3   3 MULTIPOLYGON (((35.86364 -3...
#> 4   4 POLYGON ((36.16241 -2.28064...
#> 5   5 POLYGON ((36.75994 -3.0791,...
#> 6   6 POLYGON ((35.96323 -2.77967...
#> 7   7 POLYGON ((35.2661 -1.681795...
#> 8   8 POLYGON ((39.15008 -6.77197...
#> 9   9 POLYGON ((39.05049 -6.57236...
#> 10 10 POLYGON ((39.24967 -6.87178...

# #2
species_fromr <- as.points(x = speciesr) %>% 
  st_as_sf

par(mar = c(0, 0, 0, 0))
dists_fromr %>% 
  st_geometry %>% 
  plot(col = topo.colors(n = nrow(districts)))
species_fromr %>% 
  plot(pch = 20, col = "red", add = TRUE)

Vectorizing rasters and then vectorizing back again means that you end up with lower resolution vectors if the raster is fairly coarse. You will note that this has occurred here, both in converting the rasterized districts back to polygons (#1) and the rasterized farmer counts back to points (#2). There are several things to note here. First, we used different functions for vectorizing polygons and points. Second, we piped each vectorization output to st_as_sf, because the endresult of the as.polygons/as.points functions is a SpatVector object, so the extra step coerces those to sf.

2.5.1 Questions

1. What is a primary difference between sf and s4 object classes?

2. What function should you use to read and write a multi-layer raster?

3. What are the differences between stack and brick?

4. What class of object do terra’s vectorization functions produce?

2.5.2 Code

1. Create a new raster r4, using r3 (above) as a template. Update the values of r4 using numbers randomly selected from a uniform distribution ranging between 0 and 1. Create another raster r5 by finding the values greater than 0.5 in r4. Recreate the list l with r, r2, r3, r4, and r5, and then create and plot stack s2.

2. Create a new brick b2 by applying the rast function to write s2 and to disk as b2.tif, specifying the path to your notebooks/data folder.

3. Following the steps in Chunk 8, recreate speciesr by re-rasterizing species at a 0.2 degree resolution. Plot the result.

4. Project the new 0.2 degree resolution speciesr to an Albers projection with a target resolution of 20 km (20,000 m), calling it speciesr_alb. Chunk 11 is your guide, but reproject using a bilinear rather than nearest neighbor interpolation (see ?project). Make a two panel plot comparing speciesr on the left, plotted over the unioned districts of Tanzania in grey, with speciesr_alb on the right, plotted over the unioned districts of Tanzania in grey (and transformed to Albers). Code in Chunks 8 and 11 can help.

5. Convert distsr to an sf polygons object, using as.polygons. Plot the result without coercing to st_geometry.

3.1.1 Masking

We are interested in the bioclimatic data within Tanzania’s borders, so we need to mask out the values of bioclim that fall outside of Tanzania:

# Chunk 16
test_m <- mask(x = bioclim[[1]], mask = districts)
plot(test_m, axes = FALSE, mar = c(1, 0.5, 1, 4))  

In the code above, we use the districts data to mask out the portions of bioclim (bio1) falling outside of Tanzania. Let’s apply that to the entire bioclim dataset:

# Chunk 17
bioclimz <- mask(x = bioclim, mask = districts)

set.seed(1)
ind <- sample(1:nlyr(bioclimz), size = 3)
plot(bioclimz[[ind]], axes = FALSE, nr = 1, mar = c(1, 0.5, 1, 4))

The new bioclim dataset contains all variables with the values outside of Tanzania converted to NA. The three plots are randomly selected from the layers of bioclimz, as a check to confirm that the masking was applied to all layers.

3.1.2 Cropping

If we need to chop a raster down to a smaller region, we can crop it. crop uses the extent of another object to achieve this.

# Chunk 18
bioclim1_dist72 <- crop(x = bioclimz[[1]], y = districts %>% slice(72))

plot(bioclim1_dist72, axes = FALSE, mar = c(1, 0.5, 1, 4))
plot(st_geometry(districts), add = TRUE)
districts %>% 
  st_centroid %>% 
  st_coordinates %>% 
  text(x = ., labels = 1:nrow(.))
#> Warning: st_centroid assumes attributes are constant over geometries
#> Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon = of_largest_polygon): st_centroid
#> does not give correct centroids for longitude/latitude data

The above uses the extent of district 72 (dist72) to crop the first layer of bioclimz. Note that in the plotting step we can pipe the cropped raster to plot, and also that we can label the district numbers on the plot by extracting their centroid coordinates and passing these into the “x” argument of text, and then using the nrow of the piped-in data.frame of centroid coordinates as the labels.

Cropping is important if we want to confine further analyses to a particular sub-region. However, if we just want to focus our plot to a particular region, then we can simply zoom the plot to that region (without cropping) using an extent object.

# Chunk 19
plot(bioclimz[[1]], axes = FALSE, ext = ext(districts[72, ]), 
     mar = c(1, 0.5, 1, 4))
plot(st_geometry(districts), add = TRUE)
districts %>% 
  st_centroid %>% 
  st_coordinates %>% 
  text(x = ., labels = 1:nrow(.))
#> Warning: st_centroid assumes attributes are constant over geometries
#> Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon = of_largest_polygon): st_centroid
#> does not give correct centroids for longitude/latitude data

# Chunk 20
plot_noaxes <- function(x, axes = FALSE, box = FALSE, mar = c(1, 0.5, 2, 4),
                        ...) {
  if(!class(x) %in%
     c("RasterLayer", "RasterStack", "RasterBrick", "SpatRaster")) {
    stop("This function is intended for rasters only", call. = FALSE)
  }
  par(mar = mar)
  if(class(x) %in% c("RasterLayer", "RasterStack", "RasterBrick")) {
    plot(x, axes = axes, box = axes, ...)
  } else {
    plot(x, axes = axes, mar = NA, ...)
  }
}

# Chunk 21
# bioclimz[[1]] %>% plot_noaxes  # also works in a pipeline
plot_noaxes(bioclimz[[1]])

3.1.3 Aggregating/disaggregating

To change to resolution of a raster object, you can use aggregate to make a coarser resolution, or disaggregate to decrease the resolution. To make chirpsz match the 0.1 resolution of distsr (the rasterized districts), we do this:

# Chunk 21
bioclimz1agg <- aggregate(x = bioclimz[[1]], fact = 2, fun = mean)
bioclimz1agg
#> class       : SpatRaster 
#> size        : 647, 651, 1  (nrow, ncol, nlyr)
#> resolution  : 0.01666667, 0.01666667  (x, y)
#> extent      : 29.59153, 40.44153, -11.76681, -0.9834726  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> name        :  bio01 
#> min value   : -6.075 
#> max value   : 28.650

We aggregate the first layer of bioclimz by a factor of 2 (fact = 2), taking the average of the aggregated pixels (fun = mean). Since the starting resolution was 0.05, doubling the pixel size takes it to 0.1. We could have chosen to aggregate all layers, and we could have chosen a different aggregation function (e.g. sum, max).

To disaggregate, there are two options:

# Chunk 22
bioclimz1km <- disagg(x = bioclimz[[1]], fact = 5)
#> |---------|---------|---------|---------|==============                                          
# bioclimz1km <- disagg(x = bioclimz[[1]], fact = 5, method = "bilinear") 
bioclimz1km
#> class       : SpatRaster 
#> size        : 6465, 6510, 1  (nrow, ncol, nlyr)
#> resolution  : 0.001666667, 0.001666667  (x, y)
#> extent      : 29.59153, 40.44153, -11.75847, -0.9834726  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> varname     : bioclim 
#> name        : bio01 
#> min value   : -6.65 
#> max value   : 28.65

The first option is the default one, which just breaks apart each pixel into the number of smaller new pixels specified by the factor, assigning each new pixel the same value as its larger parent. Here we specified a factor of 5, leading to an output resolution of 0.01 (~1 km), and thus 25 times the number of pixels as in the original bioclimz[[1]] (run ncell(bioclimz1km) / ncell(bioclimz[[1]]) to see).

The second way (commented out), is to set method = bilinear, which interpolates during disaggregation. In this case, disaggregate does the job using the resample function, which we will see next.

3.1.4 Resampling

Resampling is done for many reasons, but one particularly common case is as a final step following aggregation or disaggregation, when it is needed to co-register the resulting grid to another grid.

# Chunk 23
bioclim125 <- aggregate(x = bioclimz[[1]], fact = 5)  # no fun b/c default is mean
s <- c(bioclim125, speciesr)  # fails 
#> Error:
#> ! [rast] extents do not match

# Check their extent
ext(bioclim125)
#> SpatExtent : 29.59152693945, 40.46652689595, -11.77513920175, -0.983472578249996 (xmin, xmax, ymin, ymax)

ext(speciesr)
#> SpatExtent : 29.5895292580001, 40.444734823, -11.7623492139999, -0.983143016999943 (xmin, xmax, ymin, ymax)

# Very close, but not match exactly

par(mar = c(1, 1, 1, 1))
plot(ext(bioclim125), axes = FALSE, border = "blue")
plot(ext(speciesr), add = TRUE, border = "red")

After aggregating bioclimz layer 1, we try and fail to stack the result with speciesr, because (as the plot of extents show, and the warning tells) the two objects have different extents. So, another approach is needed:

# Chunk 24
speciesr_rs <- resample(x = speciesr, y = bioclim125) 
s <- c(bioclim125, speciesr_rs)  
names(s) <- c("bio1", "species_count")
plot_noaxes(s)

We resample speciesr to bioclim125, since bioclim125 has the larger extent. However, there is one more thing we should do before using these data, and that is an adjustment to speciesr. Have a look at this:

plot_noaxes(s$species_count >= 0)
plot(st_geometry(st_union(districts)), add = TRUE)
#> although coordinates are longitude/latitude, st_union assumes that they are planar

This shows that there are no places where there are 0 farmers–because the original data contributing to speciesr (the second layer in s) was derived from point locations representing the approximate locations of where species were observed. The data thus have no points where there were 0 farmers observed, thus everywhere outside of observed locations is assigned a no data value. These are not analyzed in raster operations. So we have to fix this if we want to have the full area of Tanzania in the second layer of our stack.

s$species_count[is.na(s$species_count)] <- 0
s$species_count <- mask(s$species_count, s$bio1)
plot_noaxes(s$species_count)
plot(st_geometry(st_union(districts)), add = TRUE)
#> although coordinates are longitude/latitude, st_union assumes that they are planar

Now all the areas in s$species_count where no interviews were conducted area set to 0.

We can then stack and perform subsequent analyses that draw on both layers, e.g. finding areas where annaul mean temperature exceeded 15 degree and where there was more than 1 farmer:

# Chunk 25
plot_noaxes(s$species_count > 1 & s$bio1 > 15)

# (s$bio1 > 15 & s$species_count > 1) %>% plot_noaxes  # also works

3.2.1 Global statistics

The most basic analyses to perform on rasters is to calculate global summary statistics, i.e. statistics calculated across all cell values:

# Chunk 26
global(x = bioclimz[[1]], stat = "mean", na.rm = TRUE)  # for a single variable

global(x = bioclimz[[c(5, 7, 14)]], stat = "mean", na.rm = TRUE)

summary(bioclimz[[1:3]])

global lets you calculate specific statistics from the cell values of a lone SpatRaster, or for a single, multiple, or all layers in a SpatRaster layer. summary (a generic) returns the quintile values and counts the number of NA or no data cells. Both functions by default remove NA values (which are in many if not most rasters), which is something that usually has to be specified when trying to apply these statistical functions to an ordinary vector, e.g. 

# Chunk 27
v <- values(bioclimz[[1]])
mean(v) 
#> [1] NA
mean(v, na.rm = TRUE)
#> [1] 22.43886

v is a vector of all of the values from the first layer in the bioclim brick, including its NA values. mean returns an NA if you don’t tell the function to remove NA values first. This is important to remember, because both spatial and non-spatial data often have missing values, so you will have to deal with them explicitly in many cases.

Here’s a more programmatic way of using cellStats:

# Chunk 28
# #1
cv <- function(x, na.rm) {
  cv <- sd(x, na.rm = na.rm) / mean(x, na.rm = na.rm) * 100
  return(cv)
}
bioclim_stats <- lapply(list(mean, sd, cv), function(x) {
  global(x = bioclimz, fun = eval(x), na.rm = TRUE)
})
names(bioclim_stats) <- c("mean", "sd", "cv")

# #2
bioclim_statsf <- do.call(cbind, bioclim_stats) %>% data.frame %>% 
    mutate(var = row.names(.)) %>% 
    rename(cv = global) %>% 
    mutate(type = ifelse(as.numeric(gsub("bio", "", var)) <= 11, 
                         "Temperature", "Precipitation"))

# #3
ps <- lapply(c("mean", "sd", "cv"), function(x) { 
  ggplot(bioclim_statsf, aes(x = var, y = .data[[x]], fill = type)) + 
    geom_col() + 
    facet_wrap(~type, scales = "free") + # Separate scales for temp and precip
    coord_flip() +
    scale_fill_manual(values = c("Temperature" = "#e41a1c", "Precipitation" = "#377eb8")) +
    labs(x = "Bioclim Variable", y = x, title = paste(x, "by Category")) +
    theme_minimal() +
    theme(legend.position = "none")
})
cowplot::plot_grid(plotlist = ps, ncol = 1)

In #1 we use lapply to calculate 3 different statistics over all dates in bioclimz, meaning that we get a list for each statistic as calculated over all of Tanzania.

In #2, we use a pipeline to bind the summary vectors into a data.frame (one statistic per column), and then do some gymnastics to convert the row names of the data.frame, which are constructed from the layer names of bioclimz (which are stored as “bioXX”) to a vector of variables.

We use lapply to set up four ggplot objects in #3 (note the syntax to let variable names to be passed as string, rather than unquoted), and then use cowplot::plot_grid to put those plots in a grid, using the “plotlist” argument to receive the input list of ggplots.

The results in the plots above reveal the long-term climatological structure of the landscape. The first panel shows the global mean for each layer, clearly illustrating the different scales between Temperature (BIO1–BIO11) and Precipitation (BIO12–BIO19) variables. The next two panels show two measures of spatial variability: the standard deviation (SD) and the coefficient of variation (CV). These measures summarize how much these climate factors vary across the geography of the study area.

For the temperature variables, the SD is generally low, except for BIO4 (Temperature Seasonality), which shows a distinct spike in spatial variation across the region. In the precipitation group, BIO12 (Annual Precipitation) shows a very high SD, which is expected as it tracks with the high mean values. However, the CV tells a more interesting story; while BIO12 has a high absolute variation (SD), its relative variation (CV) is actually lower than several other layers, such as BIO14 or BIO17. This indicates that while the total volume of rain varies, the most extreme spatial “patchiness” occurs during the driest periods of the year, where small local differences create significant relative spikes in variability.

Another way to summarize raster data is visually, using a histogram (terra has a generic hist function)

# Chunk 30
par(mfrow = c(1, 3))
hist(bioclimz[[1:3]], col = "blue", xlab = "mm")
#> Warning: [hist] a sample of 59% of the cells was used (of which 38% was NA)
#> Warning: [hist] a sample of 59% of the cells was used (of which 38% was NA)
#> Warning: [hist] a sample of 59% of the cells was used (of which 38% was NA)

This variant plots a histogram per layer in the SpatRaster object (but we told it to plot 1 row, 3 columns, instead of the default (2 rows, 2 columns)).

freq is another way to summarize raster values that is similar to hist but without the automatic plot.

# Chunk 31
f <- freq(bioclimz[[1]]) %>% print()
#>    layer value  count
#> 1      1    -7      3
#> 2      1    -6      6
#> 3      1    -5      4
#> 4      1    -4      3
#> 5      1    -3      7
#> 6      1    -2      7
#> 7      1    -1      9
#> 8      1     0     20
#> 9      1     1     35
#> 10     1     2     57
#> 11     1     3     58
#> 12     1     4     65
#> 13     1     5     72
#> 14     1     6     99
#> 15     1     7    103
#> 16     1     8    126
#> 17     1     9    151
#> 18     1    10    213
#> 19     1    11    529
#> 20     1    12    917
#> 21     1    13   1735
#> 22     1    14   3645
#> 23     1    15   7562
#> 24     1    16  11917
#> 25     1    17  16749
#> 26     1    18  22611
#> 27     1    19  39940
#> 28     1    20  68454
#> 29     1    21 114618
#> 30     1    22 170193
#> 31     1    23 257937
#> 32     1    24 132359
#> 33     1    25 127511
#> 34     1    26  64870
#> 35     1    27   5277
#> 36     1    28   1142
#> 37     1    29    270

Here we apply freq to a dataset with continuous values, although this function is probably best reserved for categorical rasters. However, it produces reasonable results here.

3.2.2 Local statistics

The previous section showed us how to produce statistics calculated across the entire raster. Now we will learn to calculate local, or neighborhood, statistics.

# Chunk 32
# #1
# zonemu <- zonal(x = bioclimz, z = distsr, fun = "mean")  # fails b/c extent

# #2
distsr_rs <- resample(x = distsr, y = bioclimz, method = "near")  # match extent
zonemu <- zonal(x = bioclimz, z = distsr_rs, fun = "mean", na.rm = TRUE)
head(zonemu)[, 1:5]
#>   ID    bio01     bio02     bio03    bio04
#> 1  1 18.82958 10.270139 0.6199528 183.5333
#> 2  2 20.62708 10.194445 0.6019861 202.8657
#> 3  3 21.61495  9.230065 0.6744109 145.4794
#> 4  4 21.91470 10.691425 0.6981019 126.3778
#> 5  5 19.40262 10.040972 0.6131875 183.4982
#> 6  6 21.86499 10.978851 0.6470662 179.5308

# Chunk 33
subsmat <- zonemu %>% dplyr::select(1:2)
distr_rfmu <- subst(x = distsr_rs, from = subsmat[, 1], to = subsmat[, 2])
plot_noaxes(distr_rfmu)

# Chunk 34
# #1
wmat <- matrix(1 / 9, nrow = 3, ncol = 3) 
bioclim1_focmu1 <- focal(x = bioclimz[[1]], w = wmat)

# #2
wmat <- matrix(1, nrow = 3, ncol = 3) 
bioclim1_focmu2 <- focal(x = bioclimz[[1]], w = wmat, fun = mean)

# #3
wmat <- matrix(1, nrow = 5, ncol = 5) 
bioclim1_focmu3 <- focal(x = bioclimz[[1]], w = wmat, fun = mean) 

# #4 
wmat <- matrix(1, nrow = 5, ncol = 5) 
bioclim1_focmu4 <- focal(x = bioclimz[[1]], w = wmat, fun = mean, na.rm = TRUE) 

# #5
wmat <- matrix(1 / 9, nrow = 5, ncol = 5) 
bioclim1_focmu5 <- focal(x = bioclimz[[1]], w = wmat, na.rm = TRUE) 

# plots
l <- list(bioclim1_focmu1, bioclim1_focmu3, bioclim1_focmu4, bioclim1_focmu5)
titles <- c("3X3 NAs not removed", "5X5 NAs not removed", 
            "5X5 NAs removed properly", "5X5 NAs removed improperly")

# mar: c(bottom, left, top, right) - default is c(5, 4, 4, 2) + 0.1
# oma: outer margins - default is c(0, 0, 0, 0)
par(mar = c(4, 4, 2, 1), oma = c(0, 0, 2, 0), mfrow = c(2, 2))
for(i in 1:length(l)) {
  plot_noaxes(l[[i]], main = titles[i])
}

3.2.3 Analyzing the Z dimension

If we have a SpatRaster with more than one layer, we have three dimensions (x, y, z). Often we want to analyze the values in the Z dimension (which may represent time, spectral bands, or unrelated spatial predictors in a model) without altering the x and y dimensions.

The workhorse for doing this sort of analysis is app, which allows you to apply pretty much any function to the Z-dimension of a multi-layer SpatRaster.

# Chunk 35
# #1
bioclim_zmu <- app(x = bioclimz, fun = mean)
bioclim_zsd <- app(x = bioclimz, fun = sd)
bioclim_zrng <- app(x = bioclimz, fun = range)

# #2
bioclim_zstack <- c(bioclim_zmu, bioclim_zsd, bioclim_zrng)
names(bioclim_zstack) <- c("Mean", "StDev", "Min", "Max")

plot_noaxes(bioclim_zstack)

In #1 we pass mean, sd, and range to “fun” in app. Note that range always returns two values, so app conveniently returns a two-layer SpatRaster that contains the minimum in the first layer and the maximum in the second.

In #2 we then stack the three outputs, and rename the layers something meaningful so that plot_noaxes can plot them all at once.

3.3.1 Questions

1. How do global, focal, and zonal differ from one another?

2. How do you disagg a raster with interpolation?

3. How do you run app on images with different resolutions or extents?

3.3.2 Code

1. Sometimes we will use date for analysis. Run as.Date("10-11-2017", "%m-%d-%Y"), as.Date("10-11-17", "%m-%d-%y"), as.Date("101117", "%m%d%y"), and as.Date("10112017", "%m%d%Y") to get a better sense of date vectors work. Also try out lubridate::mdy("10-11-2017") and lubridate::as_date("20171011").

2. Convert species to a 0.1 degree raster (speciesr2) that contains the count of farmers per grid cell. Use distsr as the target raster so that the extents align.

3. Use zonal on speciesr2 to calculate the total number of farmers per district (use distsr to provide the zones), and then map them back onto the districts/zones using subst. Plot the result.

4. Use focal to calculate for bioclimz[[12]] the i) standard deviation within a 3X3 and 5X5 window, and ii) the maximum value in each 3X3 and 5X5 neighborhood. Do not remove NAs. Combine the results in a multi-layer SpatRaster, as above, and then plot them using plot_noaxes.

5. Crop bioclimz[[1]] using the extent of districts[57, ], and disaggregate it to 0.01 resolution using both the default and bilinear methods. Plot the results side by side using plot_noaxes.

6. Use app with bioclimz to calculate the maximum, the coefficient of variation, and the median of all variables. Stack the results and plot with meaningful titles using plot_noaxes (hint: name the layers of the multi-layer SpatRaster), outputting plots on 1 row with 3 columns.